using System;

using Atomic.Libraries.Mathematics.LinearAlgebra.AlgLib;

namespace Atomic.Libraries.Mathematics.LinearAlgebra
{
    public sealed class CholeskyDecomposition
    {
        private CholeskyDecomposition()
        {
        }

        public static void InverseDeterminant(Matrix matrix, out Matrix inverse, out double determinant)
        {
            int n = matrix.Rows;
            if (matrix.Columns != n)
            {
                throw new ArgumentException("The matrix is not a square matrix.");
            }

            double[,] a = matrix.ToArray();

            if (!trfac.spdmatrixcholesky(ref a, n, false))
            {
                throw new ArithmeticException();
            }

            determinant = matdet.spdmatrixcholeskydet(ref a, n);

            int info = 0;
            matinv.matinvreport rep = new matinv.matinvreport();
            matinv.spdmatrixcholeskyinverse(ref a, n, false, ref info, ref rep);

            for (int i = 0; i < n; i++)
            {
                for (int j = i + 1; j < n; j++)
                {
                    a[i, j] = a[j, i];
                }
            }

            inverse = new Matrix(a);
        }
    }
}
